home *** CD-ROM | disk | FTP | other *** search
/ Die Ultimative Software-P…i Collection 1996 & 1997 / Die Ultimative Software-Pakete CD-ROM fur Atari Collection 1996 & 1997.iso / g / gnu_c / pmlsrc23.zoo / pmlsrc / atan.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-03-19  |  9.0 KB  |  357 lines
  1. /************************************************************************
  2.  *
  3.  * New version of atan function for PML library.  In general based on
  4.  * original routine by Fred Fish, but modified to use range reduction
  5.  * in a more consistent manner and to perform its task without recursive
  6.  * calls.  It does not complain when arguments are close to 0.0.
  7.  * atan() is smooth enough and its derivative there is 1.
  8.  * It also now returns PI/2 and not 0.0 when arguments are getting too big.
  9.  * Besides original function would not work when matherr() function
  10.  * would be replaced with something which always returns 1.
  11.  *
  12.  * Michal Jaegermann, December 1990
  13.  * ntomczak@ualtavm.bitnet
  14.  *
  15.  ************************************************************************/
  16.  
  17. /************************************************************************
  18.  *                                    *
  19.  *                N O T I C E                *
  20.  *                                    *
  21.  *            Copyright Abandoned, 1987, Fred Fish        *
  22.  *                                    *
  23.  *    This previously copyrighted work has been placed into the    *
  24.  *    public domain by the author (Fred Fish) and may be freely used    *
  25.  *    for any purpose, private or commercial.  I would appreciate    *
  26.  *    it, as a courtesy, if this notice is left in all copies and    *
  27.  *    derivative works.  Thank you, and enjoy...            *
  28.  *                                    *
  29.  *    The author makes no warranty of any kind with respect to this    *
  30.  *    product and explicitly disclaims any implied warranties of    *
  31.  *    merchantability or fitness for any particular purpose.        *
  32.  *                                    *
  33.  ************************************************************************
  34.  */
  35.  
  36.  
  37. /*
  38.  *  FUNCTION
  39.  *
  40.  *    atan   double precision arc tangent
  41.  *
  42.  *  KEY WORDS
  43.  *
  44.  *    atan
  45.  *    machine independent routines
  46.  *    trigonometric functions
  47.  *    math libraries
  48.  *
  49.  *  DESCRIPTION
  50.  *
  51.  *    Returns double precision arc tangent of double precision
  52.  *    floating point argument.
  53.  *
  54.  *  USAGE
  55.  *
  56.  *    double atan (x)
  57.  *    double x;
  58.  *
  59.  *  REFERENCES
  60.  *
  61.  *    Fortran 77 user's guide, Digital Equipment Corp. pp B-3
  62.  *
  63.  *    Computer Approximations, J.F. Hart et al, John Wiley & Sons,
  64.  *    1968, pp. 120-130.
  65.  *
  66.  *  RESTRICTIONS
  67.  *
  68.  *    The maximum relative error for the approximating polynomial
  69.  *    is 10**(-16.82).  However, this assumes exact arithmetic
  70.  *    in the polynomial evaluation.  Additional rounding and
  71.  *    truncation errors may occur as the argument is reduced
  72.  *    to the range over which the polynomial approximation
  73.  *    is valid, and as the polynomial is evaluated using
  74.  *    finite-precision arithmetic.
  75.  *    
  76.  *  PROGRAMMER
  77.  *
  78.  *    Fred Fish
  79.  *
  80.  *  INTERNALS
  81.  *
  82.  *    Computes arctangent(x) from:
  83.  *
  84.  *        (1)    If x < 0 then negate x, perform steps
  85.  *            2, 3, and 4, and negate the returned
  86.  *            result.  This makes use of the identity
  87.  *            atan(-x) = -atan(x).
  88.  *
  89.  *        (2)    If argument x > 1 at this point then
  90.  *            test to be sure that x can be inverted
  91.  *            without underflowing.  If not, reduce
  92.  *            x to largest possible number that can
  93.  *            be inverted, issue warning, and continue.
  94.  *            Perform steps 3 and 4 with arg = 1/x
  95.  *            and subtract returned result from
  96.  *            pi/2.  This makes use of the identity
  97.  *            atan(x) = pi/2 - atan(1/x) for x>0.
  98.  *
  99.  *              (2a)
  100.  *                      Modification - numbers in range [1, tan((5/12)pi)]
  101.  *                      are not iverted as above but a range is reduced
  102.  *                      to [-tan(pi/12),tan(pi/12)] with result modified
  103.  *                      correspondingly ( -- mj)
  104.  *
  105.  *        (3)    At this point 0 <= x <= 1.  If
  106.  *            x > tan(pi/12) then perform step 4
  107.  *            with arg = (x*sqrt(3)-1)/(sqrt(3)+x)
  108.  *            and add pi/6 to returned result.
  109.  *                      (modification - same formula, different rendition --mj)
  110.  *            This last transformation maps arguments
  111.  *            greater than tan(pi/12) to arguments
  112.  *            in range 0...tan(pi/12).
  113.  *
  114.  *        (4)    At this point the argument has been
  115.  *            transformed so that it lies in the
  116.  *            range -tan(pi/12)...tan(pi/12).
  117.  *                      (The check below eliminated - it really
  118.  *                      doesn't make sense.  Approximation is
  119.  *                      surprisingly robust anyway -- mj)
  120.  *            Since very small arguments may cause
  121.  *            underflow in the polynomial evaluation,
  122.  *            a final check is performed.  If the
  123.  *            argument is less than the underflow
  124.  *            bound, the function returns x, since
  125.  *            atan(x) approaches asin(x) which
  126.  *            approaches x, as x goes to zero.
  127.  *
  128.  *        (5)    atan(x) is now computed by one of the
  129.  *            approximations given in the cited
  130.  *            reference (Hart).  That is:
  131.  *
  132.  *            atan(x) = x * SUM [ C[i] * x**(2*i) ]
  133.  *            over i = {0,1,...8}.
  134.  *
  135.  *            Where:
  136.  *
  137.  *            C[0] =    .9999999999999999849899
  138.  *            C[1] =    -.333333333333299308717
  139.  *            C[2] =    .1999999999872944792
  140.  *            C[3] =    -.142857141028255452
  141.  *            C[4] =    .11111097898051048
  142.  *            C[5] =    -.0909037114191074
  143.  *            C[6] =    .0767936869066
  144.  *            C[7] =    -.06483193510303
  145.  *            C[8] =    .0443895157187
  146.  *            (coefficients from HART table #4945 pg 267)
  147.  *
  148.  */
  149.  
  150. #include <stdio.h>
  151. #include <math.h>
  152. #include "pml.h"
  153.  
  154. #if !defined (__M68881__) && !defined (sfp004)    /* mjr++        */
  155.  
  156. static char funcname[] = "atan";
  157.  
  158. static double atan_coeffs[] = {
  159. /*  .9999999999999999849899, */        /* P0 must be first    */
  160.     .99999999999999998,                 /* watch out - bug in gcc 1.37
  161.                                            at least for ST */
  162.                                         /* longer representations will burn
  163.                        you badly -- mj */
  164.     -.333333333333299308717,
  165.     .1999999999872944792,
  166.     -.142857141028255452,
  167.     .11111097898051048,
  168.     -.0909037114191074,
  169.     .0767936869066,
  170.     -.06483193510303,
  171.     .0443895157187                /* Pn must be last    */
  172. };
  173.  
  174. #define LAST_BOUND 0.2679491924311227074725    /*  tan (PI/12)        */
  175. #define MID_BOUND  1.0                /*  tan (3 * PI/12)    */
  176. #define HI_BOUND   3.7320508075688772935274    /*  tan (5 * PI/12)    */
  177.  
  178. #define INV_ROOT3  0.5773502691896257645
  179. #define THIRDPI    1.0471975511968977461
  180.  
  181.     struct exception xcpt;
  182.  
  183. double atan (x)
  184. double x;
  185. {
  186.     short neg_flag = 0;
  187.  
  188.     xcpt.retval = 0.0;
  189.     if (x < 0.0) {
  190.     x = -x;
  191.     neg_flag |= 1;
  192.     }
  193.     if (x > MID_BOUND) {
  194.     if (x > HI_BOUND) {
  195.         xcpt.retval = HALFPI;
  196.         if (x >= MAXDOUBLE) {
  197.         if (!matherr (&xcpt)) {
  198.             xcpt.type = UNDERFLOW;
  199.             xcpt.name = funcname;
  200.             xcpt.arg1 = x;
  201.             fprintf (stderr, "%s: UNDERFLOW error\n", funcname);
  202.             errno = EDOM;
  203.         }
  204.         x = 0.0;
  205.         }
  206.         else {
  207.         x = -1.0 / x;
  208.         }
  209.     }
  210.     else {  /* MID_BOUND <= x <= HI_BOUND */
  211.         xcpt.retval = THIRDPI;
  212.         x = INV_ROOT3 -  4.0 / (SQRT3 + x + x + x);
  213.     }
  214.     }
  215.     else { /* (x <= MID_BOUND) */
  216.     if (x > LAST_BOUND) {
  217.         xcpt.retval = SIXTHPI;
  218.         x = SQRT3 - 4.0 / (SQRT3 + x);
  219.     }
  220.     }
  221.     /* range reduced to [-LAST_BOUND, LAST_BOUND] */
  222.     if (x != 0.0) { /* polynomial approximation */
  223.                 /* underflow in poly is not of a big concern here */
  224.     x *= poly (((int)sizeof (atan_coeffs) / (int)sizeof(double)) - 1,
  225.             atan_coeffs, x * x);
  226.     xcpt.retval += x;
  227.     }
  228.     if (neg_flag)
  229.         xcpt.retval = -xcpt.retval;
  230.     return (xcpt.retval);
  231. }
  232. #endif    /* !__M68881__ || !sfp004    */
  233. #ifdef    sfp004
  234.  
  235. __asm("
  236. comm =     -6
  237. resp =    -16
  238. zahl =      0
  239. ");    /* end asm    */
  240.  
  241. #endif    sfp004
  242. #if defined (__M68881__) || defined (sfp004)
  243.     __asm(".text; .even");
  244.  
  245. # ifdef    ERROR_CHECK
  246.  
  247.     __asm("
  248.  
  249. _Overflow:
  250.     .ascii \"OVERFLOW\\0\"
  251. _Domain:
  252.     .ascii \"DOMAIN\\0\"
  253. _Error_String:
  254.     .ascii \"atan: %s error\\n\\0\"
  255. .even
  256.  
  257.  
  258. | pml compatible atangent
  259. | m.ritzert 7.12.1991
  260. | ritzert@dfg.dbp.de
  261. |
  262. |    /* NAN  = {7fffffff,ffffffff}        */
  263. |    /* +Inf = {7ff00000,00000000}        */
  264. |    /* -Inf = {fff00000,00000000}        */
  265. |    /* MAX_D= {7fee42d1,30773b76}        */
  266. |    /* MIN_D= {ffee42d1,30773b76}        */
  267.  
  268. .even
  269. double_max:
  270.     .long    0x7fee42d1
  271.     .long    0x30273b76
  272. double_min:
  273.     .long    0xffee42d1
  274.     .long    0x30273b76
  275. NaN:
  276.     .long    0x7fffffff
  277.     .long    0xffffffff
  278. p_Inf:
  279.     .long    0x7ff00000
  280.     .long    0x00000000
  281. m_Inf:
  282.     .long    0xfff00000
  283.     .long    0x00000000
  284. ");
  285. # endif    ERROR_CHECK
  286.  
  287. __asm("
  288. .even
  289.     .globl _atan
  290. _atan:
  291.     ");    /* end asm    */
  292.  
  293. #endif    /* __M68881__ || sfp004    */
  294. #ifdef    __M68881__
  295.  
  296.     __asm("
  297.     fatand    a7@(4), fp0    | atan
  298.     fmoved    fp0,a7@-    | push result
  299.     moveml    a7@+,d0-d1    | return_value
  300.     ");    /* end asm    */
  301.  
  302. #endif    __M68881__
  303. #ifdef    sfp004
  304.     __asm("
  305.     lea    0xfffa50,a0
  306.     movew    #0x540a,a0@(comm)    | specify function
  307.     cmpiw    #0x8900,a0@(resp)    | check
  308.     movel    a7@(4),a0@        | load arg_hi
  309.     movel    a7@(8),a0@        | load arg_low
  310.     movew    #0x7400,a0@(comm)    | result to d0
  311.     .long    0x0c688900, 0xfff067f8    | wait
  312.     movel    a0@,d0
  313.     movel    a0@,d1
  314.     ");    /* end asm    */
  315.  
  316. #endif    sfp004
  317. #if defined (__M68881__) || defined (sfp004)
  318. # ifdef    ERROR_CHECK
  319.     __asm("
  320.     lea    double_max,a0    |
  321.     swap    d0        | exponent into lower word
  322.     cmpw    a0@(16),d0    | == NaN ?
  323.     beq    error_nan    |
  324.     swap    d0        | result ok,
  325.     rts            | restore d0
  326. ");
  327. #ifndef    __MSHORT__
  328. __asm("
  329. error_nan:
  330.     moveml    a0@(24),d0-d1    | result = +inf
  331.     moveml    d0-d1,a7@-
  332.     movel    #62,_errno    | NAN => errno = EDOM
  333.     pea    _Domain        | for printf
  334. ");
  335. #else    __MSHORT__
  336. __asm("
  337. error_nan:
  338.     moveml    a0@(24),d0-d1    | result = +inf
  339.     moveml    d0-d1,a7@-
  340.     movew    #62,_errno    | NAN => errno = EDOM
  341.     pea    _Domain        | for printf
  342. ");
  343. #endif    __MSHORT__
  344. __asm("
  345. error_exit:
  346.     pea    _Error_String    |
  347.     pea    __iob+52    |
  348.     jbsr    _fprintf    |
  349.     addl    #12,a7        |
  350.     moveml    a7@+,d0-d1
  351.     ");
  352. # endif    ERROR_CHECK
  353.  
  354. __asm("rts");
  355.  
  356. #endif /* __M68881__ || sfp004    */
  357.